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A large step in the validation of Computational Fluid Dynamics (CFD) for Supersonic 
Retropropulsion (SRP) is shown through the comparison of three Navier-Stokes solvers 
(DPLR, FUN3D, and OVERFLOW) and wind tunnel test results. The test was designed 
specifically for CFD validation and was conducted in the Langley supersonic 4’x4’ Unitary 
Plan Wind Tunnel and includes variations in the number of nozzles, Mach and Reynolds 
numbers, thrust coefficient, and angles of orientation. Code-to-code and code-to-test com- 
parisons are encouraging and possible error sources are discussed. 


Nomenclature 


a 

Angle of Attack 

Cp 

Thrust Coefficient, T/qA re f 

0 

Roll Angle 

L 

Model Length, 10.2 in 

J^-ref 

Reference Area, 19.63 in 2 

L r ef 

Reference Length, 5.0 in 

C A, aero 

, Aerodynamic Contribution to Axial Force 

Q 

Dynamic Pressure, psi 


Coefficient 

R 

Model Radius, 2.5 in 

C A, total 

r Total Axial Force Coefficient, Cp + C a, aero 

r 

Radial Coordinate 

Cm 

Pitching Moment Coefficient 

T 

Thrust, Ibf 

Cn 

Normal Force Coefficient 

X 

Axial Coordinate 

Cp 

Pressure Coefficient 




I. Introduction 

S upersonic Retropropulsion (SRP) is a viable means for deceleration of high mass vehicles entering into 
the Martian atmosphere. 1, 2 ,3, 4,5,6 Previous methods of deceleration are not scalable for exploration type 
vehicles which can potentially weigh tens of metric tons. Since ground and flight testing of SRP at entry 
conditions can be difficult and cost-prohibitive, the development of this enabling technology can be enhanced 
with the ability to predict the flowfield numerically using Computational Fluid Dynamics (CFD). 

SRP results in a complex flow structure involving shocks, shear layers, recirculation and stagnation 
regions, which makes validation of the CFD methods a high priority. The validation process includes using 
multiple CFD codes to compare to historical 7 and recent wind tunnel tests. 8 Through code-to-code and 
code-to-test comparisons, best practices in gridding, numerical method selection, and solution advancement 
are established, and validity is added to the CFD methods. With continuing validation, confidence can be 
built for using CFD to predict Mars entry conditions. 9 
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Three CFD codes are being applied to SRP: DPLR (Data Parallel Line Relaxation), 10 FUN3D (Fully 
Unstructured Navier-Stokes Three-Dimensional), 11,12 and OVERFLOW (OVERset grid FLOW solver). 13 
The codes all solve the Reynolds Averaged Navier-Stokes equations, but differ in implementation, grid type, 
and numerical methods. The focus of this paper will be on the comparison of the CFD codes to a recent 
wind tunnel test which was designed primarily for CFD validation. The experiment was conducted by the 
NASA Exploration Technology Development Program in the Langley supersonic 4’x4’ Unitary Plan Wind 
Tunnel in July, 2010. 14, 15,16 The cases that will be presented all have a free stream Mach number of 4.6 and 
a Reynolds number of 1.5E+06 per foot, but vary by the number of nozzles (0, 1, 3, or 4 nozzles), thrust 
coefficient (Ct = T/qA re f = 2,3), angle of attack (0, 12, 16, and 20 degrees), and roll angle (0 and 180 
degrees). Time-accurate CFD simulations were conducted due to the inherent unsteadiness of the flowfields. 

Qualitative comparisons of the flow structure will be made by comparing CFD to high-speed Schlieren 
images, and quantitative comparisons will be made by comparing averaged surface pressure with pressure 
tap data from the tunnel. Unsteady shedding frequencies of the CFD solutions are also compared to high- 
frequency pressure gauges from the test. 

This paper will first introduce SRP, the CFD codes, and the wind tunnel test. Then it will present 
code-to-code and code-to-test comparisons, discuss the results including modeling strengths and weaknesses, 
and offer conclusions of the study. 


II. SRP Flow Structure 


Figure 1 shows a representative wind tunnel model that em- 
ploys SRP along with the resulting flowfield. The vehicle is a 
sphere cone with a single nozzle in the center and a sting for 
mounting it in the tunnel. In supersonic flow a bow shock forms. 
As an opposing jet is initiated, the bow shock is pushed up- 
stream as the apparent body frontal area is increased with the 
appearance of a barrel plume. The barrel plume contains free 
shear layers as well as a terminal shock. Between the terminal 
and bow shocks is an interface or contact surface where opposing 
velocities meet and go to zero. With the barrel plume, recircu- 
lation regions appear as well as a triple-point, where it can be 
said that three types of flow meet- supersonic jet flow, subsonic 
shock-layer flow, and subsonic recirculating flow. 

This is an example of a single jet flow at a relatively high 
thrust coefficient. Other modes exist depending on jet-to-free 
stream pressure ratios, including a long-jet penetration mode at 
low ratios. With multiple nozzles, interaction between barrel 
plumes may also exist, and the apparent body frontal area can 
change depending on the location of the nozzle on the model face, 
and orientation angle of the nozzles. 17,18 

III. Description of CFD Solvers 

The three solvers applied to the SRP problem differ in im- 
plementation, grid type, and numerical methods. DPLR and 
FUN3D are finite- volume while OVERFLOW is finite-difference. 
DPLR uses cell-centered structured overset grids, while OVER- 
FLOW uses node-centered structured overset grids. FUN3D em- 
ploys node centered unstructured grids. With these differences 
between codes, much is to be learned through code-to-code com- 
parison when applying them to a single set of problems. 

This paper builds from reference [9] , which focuses on a single 
run of the Langley Unitary Wind Tunnel test. In reference [9], 
an in-depth look at grid and time resolution, observed order of 
accuracy, and the establishment of best practices is described in 


Supersonic 

Freestream 





Figure 1. SRP flow structure description 
diagrams. 
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detail. The following will be brief descriptions of each code including the numerical methods used by each 
in this study. 

III. A. DPLR 

The DPLR CFD code is a parallel, structured multi-block, finite-volume code with overset grid capability 
that solves the Reynolds- Averaged Navier- Stokes (RANS) equations for continuum flow, including finite-rate 
chemistry and thermal non-equilibrium. In the present study, the thermally- and calorically-perfect RANS 
equations for air are solved implicitly with lst-order time accuracy. Inviscid fluxes are formed via a modified 
Steger- Warming flux vector splitting 19 with 3rd-order Monotone Upwind Scheme for Conservation Laws 
(MUSCL) extrapolation 20 subject to a minmod limiter 21 and 2nd-order trapezoidal rule flux integration. 
The viscous terms are computed with 2nd-order spatial accuracy with a central difference approach. For 
the present analysis, the Menter’s Shear-Stress Transport (SST) turbulence model 22 was employed with a 
vorticity-based production term and no compressibility corrections. 

III.B. FUN3D 

FUN3D contains a node-based finite-volume unstructured flow solver. For this study, FUN3D was run 
with a selectively-dissipative version of the Low-Dissipation Flux Splitting Scheme (LDFSS) inviscid flux 
function 23 and a modified Van Albada 24 limiter according to White. 25 Direct Eddy Simulation (DES) tur- 
bulence modeling was used with Spalart-Allmaras (SA) 26 as the submodel. Cases were also completed with 
the compressible RANS equations loosely coupled to Menter’s SST turbulence model with a vorticity-based 
production term. Node-based conservative variables are computed by driving a 2nd-order accurate spatial 
residual to stead-state with a point-implicit iterative method. A modified, optimum 2nd-order backward 
difference formula (BDF) scheme is used in conjunction with a temporal error controller that assured design 
order. 27 

III.C. OVERFLOW 

OVERFLOW 2 is an implicit Reynolds Averaged Navier-Stokes (RANS) flow solver that utilizes structured 
overset grids. For the current work, the HLLE++ numerical flux function 28 with the Van Albada limiter 
was used for spatial terms, and the Symmetric Successive Over Relaxation (SSOR) algorithm with dual 
time stepping using Newton sub-iterations for temporal terms. All viscous terms were included, and the 
SST turbulence model with strain-based production was employed with Wilcox’s realizability constraint. 
The overall scheme is 2nd-order accurate in space and time. The calculation of the inviscid fluxes for both 
the flow solver and the turbulence model use 3rd-order accurate MUSCL reconstruction and 2nd-order flux 
quadrature. For the single and triple nozzle configurations, DES turbulence modeling with SST as the 
submodel was used, and for the baseline and quad nozzle configurations, RANS SST turbulence modeling 
was used. 


IV. SRP Wind Tunnel Test 

Test 1853 in the Langley supersonic Unitary Plan Wind Tunnel (LaUPWT) was designed specifically 
for SRP CFD validation. The model was a 70° sphere-cone forebody with a cylindrical side body with a 
5 inch diameter. The model included four nozzles which could be plugged to offer a 0, 1, 3, or 4 nozzle 
configuration. One nozzle was located at the center of the forebody, and three were oriented at a constant 
1/2 radial coordinate every 120°. Air was used as both the freestream and jet gases. The test data included 
high speed Schlieren movies and pressure readings from 172 originally planned taps including nine 40 kHz 
pressure transducers. From these data, it is possible for qualitative comparisons with CFD of flow structure 
and unsteady behavior with the Schlieren movies, as well as averaged surface pressure comparisons with the 
pressure taps. The high frequency pressure transducers also allows frequency comparisons to measure the 
accuracy of the CFD unsteadiness. 

The test run matrix included four nozzle configurations: 0, 1 (center), 3 (peripheral), and 4 (center and 
peripheral). Three Mach numbers were tested: 2.4, 3.5, and 4.6. The Re/ft for the two lower Mach numbers 
was 1M, and for Mach 4.6, Re/ft was 1.5M. Thrust coefficients ranged from 0.5 to 3, with isolated runs up 
to 6. Angles of attack swept from -8° to 20°. A full description of the test can be found in reference [14]. 
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Figure 2 is a snapshot of the model in the triple nozzle configuration in the test section. Figure 3 shows 
diagrams of the model face and side-body pressure tap layouts. 



Figure 2. SRP test model in the 3 nozzle configuration in the NASA Langley UPWT 4x4 supersonic tunnel (test 
section 2). 



Figure 3. Pressure tap layout on test model. Left image is the front face of the model and the right image is the side 
body. Filled circles represent the high frequency pressure transducers. 


IV. A. Wind Tunnel Uncertainty 

Wind tunnel uncertainty was calculated using an entirely statistical approach based on comparing multiple 
measurements of surface pressure coefficient in a specially designed sampling technique. 29 The uncertainty 
contains contributions from random sources (instrumentation drift, hysteresis, etc), flowfield non- uniformity, 
and model geometry /instrumentation uncertainties. The final pressure coefficient uncertainty for Mach 4.6 
is seen in Table 1, where it is seen that flowfield non- uniformity is the largest contributor. For the CFD 
comparisons, 3cr uncertainty bars were applied to the test data. 

The uncertainty analysis was conducted on the baseline (no jet) configuration. As such, it does not 
include contributions introduced by the jets or the unsteadiness in the plumes. Although these contributions 
cannot be quantified directly, scatter in the results of repeat blowing cases were largely contained within 
certainty limits. 29 
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Table 1. Wind Tunnel Uncertainty 


Mach 

Source 

<j 

% Total 

4.6 

Random 

0.00294 

15 

4.6 

Flow Field 

0.00637 

71 

4.6 

Geometry 

0.00276 

13 

4.6 

Total 

0.00754 

100 


V. Results 

A total of 18 cases were simulated with the CFD solvers. They focused on the higher Mach number 
(4.6) to decrease the probability of wall interference in the tunnel which had large boundary layers in the 
test section. The cases vary in angle of attack, roll angle, thrust coefficient, and the number of nozzles (see 
Table 2). 

Images of the grids used for each code are seen in Figure 4. Each grid has a different topology but all 
focus refinement in the shock interaction region. 

Comparisons of results between the three CFD flow solvers and the wind tunnel data are shown in 
Figures 5 through 26. A general explanation of the plots will be offered here and a full discussion of the 
results will be presented later. 

The black and white images labeled as “TEST” are instantaneous images gathered from the high speed 
Schlieren system from the LaUPWT test. Colored flow visualizations from CFD codes show either Mach or 
Cp contours. Black and white flow visualizations from the CFD codes are simulated Schlieren or shadowgraph 
representations generated through volume integration as described by Yates. 30 The exception would be the 
OVERFLOW run 165 cases which are symmetry plane slices of the log of the density gradient magnitude 
colored in grayscale. To account for </> = 180° rotations, negative angles of attack were used for the FUN3D 
results. 

Columns in the flow visualizations montages are snapshots in time (labeled Time 1, Time 2, Time 3, and 
Time 4), but they do not correlate between the CFD simulations or test Schlieren. Each point in time was 
chosen to best represent the behavior observed in the flowfield. 

In the line plots, comparisons between the codes and experimental values of the averaged surface Cp 
are shown using two coordinate systems. The first shows pressures on the model face in radial coordinates 
normalized by the model radius. Data from both 0 = 0° and </> = 180° (see Figure 3) are shown on the 
same plot. The second set of plots show pressures down the length of the body. Here the X coordinates are 
normalized by the model length. 

Forces and moments as a function of time are also plotted as a code-to-code comparison. The test did 
not include a balance and integration of the pressure tap values only yield averaged aerodynamic loads. 
Total axial force is plotted, which is a summation of the aerodynamic and thrust components to axial force. 
Through these plots it is possible to compare the level of unsteadiness between the codes as well as gain an 
understanding of how much the unsteadiness may influence the stability of a vehicle at these conditions. The 
force and moment integration was conducted on the model face and side; it did not include the model base. 
The reference length and area used in nondimensionalizing the coefficients are seen in Table 3; the reference 
length is the diameter of the model, and the reference area is the projected frontal area. The moment 
reference center was the nose of the model. The time period shown in the force and moment coefficient plots 
do not represent the total run time of the CFD codes but was based on the code which was run for the 
smallest amount of time. 
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Table 2. CFD Run Matrix 


Run Number 

Nozzles 

Ct 

Angles of Attack (<a) 

Roll Angle (0) 

283 

0 

N/A 

0°, 12°, 20° 

0° 

165 

1 

2 

0°, 12°, 20° 

o 

o 

00 
T— 1 

262 

3 

3 

0°, 12°, 16° 

0° 

263 

3 

3 

0°, 12°, 16° 

o 

o 

00 
T— 1 

307 

4 

2 

0°, 12°, 20° 

0° 

311 

4 

2 

0°, 12°, 20° 

o 

o 

00 
T— 1 


Table 3. Reference Quantities for Force and Moment Coefficients 

Reference Length ( L re f ) Reference Area ( A re f ) 

5.0 in 19.63 in 2 



Figure 4. Symmetry plane slices of the grids used for each code. 
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V.A. Run 283- 0 Nozzles 


Run 283 is a baseline zero nozzle case resulting in a steady supersonic blunt-body flow, which is simpler than 
the powered cases. As expected, each code compared well to each other and to the test in shock standoff 
distances and surface Cp. The DPLR visualizations show very little flowfield upstream of the bow shock 
which demonstrates the grid adaptation feature in the DPLR code. 

Overall, the CFD codes predict higher pressure on the model face than the experiment, most of which is 
contained inside the error bars. However, the pressure at the nose is an area of deviation between the CFD 
codes and the test data for a = 12° and a = 20°. The Rayleigh pitot formula gives the maximum Cp for 
Mach 4.6 flow to be 1.803. This correlates to the Cp at the nose for a = 0°. The test data and OVERFLOW 
results show a slightly lower value than 1.803 while the DPLR and FUN3D results are closer to the theoretical 
value. This may show that the test data is biased slightly low, perhaps due to flow non-uniformity in the 
tunnel, but this deviation at the nose for a = 0° is still within the uncertainty bars. 



Figure 5. Flow visualizations of run 283, baseline zero nozzles. 
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Figure 6. Pressure coefficient on model surface for run 283, baseline zero nozzles. 
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V.B. Run 165, Single Nozzle, Cp 2 


Run 165 o = 0° was used to establish the best practices for each code. 9 It was chosen because of its periodic 
unsteadiness, which simplified comparisons. The unsteadiness can be described as an oscillation of the triple 
point (described three dimensionally as an annular ring) which created pressure waves that propagate up to 
the bow shock causing a minimal effect, and also propagate to the model. The waves reflect off the model 
face and in turn off the barrel plume shear layers which causes another oscillation of the triple point. 

Each code captured this unsteadiness to a different degree, the differences in turbulence modeling and 
grid refinement being the key contributors to the deviations. The level of unsteadiness can be seen in the 
force and moment plots in Figure 10, OVERFLOW and FUN3D having larger oscillation amplitudes than 
DPLR, with the DPLR case trending toward a steady solution. 

An effect of capturing the described unsteady effects can be noted in the model face surface Cp. The 
codes that captured the unsteadiness captured the pressure wave as it reflected off the model face, and in 
turn a higher averaged pressure near the nozzle was predicted. This trend is also seen in the wind tunnel 
data. However, the tunnel uncertainties are large enough to envelope the CFD data that did not capture 
the unsteady effects. 

For a = 12°, the flow unsteadiness becomes less periodic. A shedding from the windward side of the triple 
point occurred which was again captured to different degrees between the codes. More deviation between the 
codes and test is seen in the Cp plots. Ah codes predict pressures on the side body within test uncertainties, 
but only OVERFLOW does for the complete model face. More deviation on the windward side of the face 
is seen, especially as r/R approaches 1, or approaching the shoulder. FUN3D overpredicts pressure near the 
shoulder while DPLR underpredicts. 

For a = 20° , the how unsteadiness loses ah periodicity. Large shedding occurs randomly, and a large run 
time was required to obtain meaningful averages from the CFD cases (the same can be said for ah randomly 
unsteady cases in this paper) . The CFD results mostly fall inside the tunnel uncertainty for averaged surface 
Cp even though different behaviors are seen in the simulations. The DPLR results appear to have a periodic 
unsteadiness in the force and moment plots while the OVERFLOW and FUN3D results show a more chaotic 
behavior. 



Figure 7. Flow visualizations of run 165, ol 0°, single nozzle, Ct 2. 
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OVERFLOW FUN3D DPLR TEST OVERFLOW FUN3D DPLR TEST 



Figure 8. Flow visualizations of run 165, ol 12°, single nozzle, Ct 2. 


Time 1 


Time 2 


Time 3 


Time 4 


Time 1 Time 2 Time 3 Time 4 



Figure 9. Flow visualizations of run 165, a: 20°, single nozzle, Ct 2. 
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(a) Model surface Cp 



(b) Model C a, total, C N , and C m 

Figure 10. Comparison of averaged CFD and test surface pressure coefficient as well as code-to-code comparisons of 
C a, total, C]\r , and Cm as a function of time. Run 165, single nozzle, Ct 2. 
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V.C. Run 262, Triple Nozzle, Ct 3, (p O' 


Run 262 showed random unsteadiness in the plume to freestream interaction. The unsteadiness was shown 
in the bow shock behavior, which oscillated to different degrees throughout the testing time window. Qual- 
itatively, the FUN3D and OVERFLOW results most resemble the test Schlieren; however, that did not 
guarantee a good comparison in surface pressure. In the DPLR simulated shadowgraph images, the skew- 
ness in the area just aft of the model shoulder is due to grid effects of an overset region and a drop in 
resolution (see Figure 4). 

The level of unsteadiness for each code is apparent in the force and moment plots seen in Figure 14. 
Large deviation is seen in the total axial force between the codes. The thrust component of the coefficient 
overpowers the aerodynamic effects, and as such, small deviations in properly simulating the thrust coefficient 
become apparent. The computed Ct from the test was 2.937, which is closely simulated by FUN3D (2.99) 
and OVERFLOW (2.975), DPLR simulated a larger Ct (3.109). 

For the surface pressure, large deviations between the codes are seen at the model nose. None of the 
codes properly predict the pressure at the nose for a ® 0°, OVERFLOW comes close for a = 12°, and better 
agreement is seen for a = 16°. Differences in simulated jet expansion and jet-to-jet interactions may be the 
source of these deviations. 

For both a = 0° and a = 12°, FUN3D overpredicts the pressure for r/R < 0.4 and for the entire model 
face for a = 12°. The rise in pressure near the nozzle at r/R ~ 0.35 is captured in behavior and value by 
OVERFLOW, and by behavior for FUN3D. Large differences are seen on the windward side of the model 
face near the shoulder at non-zero angles of attack. In that region, DPLR and FUN3D properly predict 
the pressure for a = 16° with OVERFLOW overpredicting. At a = 12° a different trend is seen with 
OVERFLOW giving the best prediction of the three codes. On the model side-body, DPLR is consistently 
higher in average pressure than FUN3D and OVERFLOW with varying levels of agreeance with tunnel data 
for all codes. 


Time 1 Time 2 Time 3 Time 4 



Figure 11. Flow visualizations of run 262, ol 0°, <fi 0°, triple nozzle, Ct 3. 
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OVERFLOW FUN3D DPLR TEST OVERFLOW FUN3D DPLR TEST 


Time 1 Time 2 Time 3 Time 4 



Figure 12. Flow visualizations of run 262, ol 12°, <fi 0°, triple nozzle, Ct 3. 


Time 1 Time 2 Time 3 Time 4 



Figure 13. Flow visualizations of run 262, ol 16°, <fi 0°, triple nozzle, Ct 3. 
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prppi. 




(a) Model surface Cp 



(b) Model C A , total, C N: and C m 

Figure 14. Comparison of averaged CFD and test surface pressure coefficient as well as code-to-code comparisons of 
C a, total, Cn , and Cm as a function of time. Run 262, triple nozzle, Ct 3, <fi 0°. 
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V.D. Run 263, Triple Nozzle, Ct 3, (p 180 


Run 263 differs from run 262 only in roll angle, run 263 having the nozzle which lays in the symmetry plane 
facing windward for non-zero angles of attack. Little difference from run 262 was seen or expected at a = 0° . 
However, for a = 12° and a = 16°, run 263 displayed steadier behavior than run 262. The windward side 
of the flowfield demonstrated a firmly placed bow shock, and little differences were seen in the barrel plume 
or termination shocks. Shedding from the jets did occur on the leeward side which pushed the bow shock 
in that area back in a semi-oscillatory manner. The shedding was then pushed downstream, which for the 
most part did not effect the model surface. 

The steadier flow behavior for non-zero angles of attack was captured the best by FUN3D. DPLR and 
OVERFLOW showed unsteady behaviors similar to run 262. 

Pressure at the nose is again a weak point in the CFD simulations. Large deviation is seen between the 
codes and none of the codes properly match the tunnel data. The pressure on the model face for a = 0° 
was overpredicted by FUN3D. For the non-zero angles of attack on the model side shell, the DPLR results 
compare the best with tunnel data on the windward side while FUN3D and OVERFLOW underpredict for 
a = 12°, and OVERFLOW underpredicts for a = 16°. 
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Figure 15. Flow visualizations of run 263, ct 0°, <fi 180°, triple nozzle, Ct 3. 



Time 1 Time 2 Time 3 Time 4 
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OVERFLOW FUN3D DPLR TEST OVERFLOW FUN3D DPLR TEST 



Figure 16. Flow visualizations of run 263, ex 12°, ej) 180°, triple nozzle, Ct 3. 


Time 1 Time 2 Time 3 Time 4 



Figure 17. Flow visualizations of run 263, ex 16°, ef> 180°, triple nozzle, Ct 3. 
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(a) Model surface Cp 




(b) Model C A , total, C N: and C m 

Figure 18. Comparison of averaged CFD and test surface pressure coefficient as well as code-to-code comparisons of 
C a, total, Cn , and Cm as a function of time. Run 263, triple nozzle, Ct 3, <fi 180°. 
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V.E. Run 307, Quad Nozzle, Cp 2, (p O' 


Run 307 at a = 0° demonstrated a steadier mode than seen in any of the triple nozzle cases. The flow was 
similar to run 165 a = 0° in that periodic oscillations of the triple points occurred and the bow shock is in a 
steadier state relatively close to the body. However, the flow physics are more complex now with the added 
plume to plume interactions. For a = 12° and a = 20°, similar chaotic unsteadiness to run 262 is seen with 
fluctuations in the bow shock. 

For a = 0°, FUN3D and OVERFLOW properly simulate the steadier flowfield with the short bow 
shock standoff distance. The DPLR simulation shows behavior of a larger bow shock standoff distance as it 
approaches a steady state. The bow shock standoff distance does not seem to influence the surface pressure, 
however, since the plume structure is large enough to block oncoming freestream flow from the model. The 
pressure on the model is low, and all codes fall within experimental uncertainties with the exception of 
OVERFLOW below the center nozzle. 

For a = 12° and a = 20°, all codes predict a higher level of unsteadiness, FUN3D picking up the most 
fluctuations For the most part, all codes predict Cp well on the side of the model, except OVERFLOW 
underpredicts aft of the shoulder for a = 12°, and FUN3D underpredicts aft of the shoulder for a = 20°. On 
the model face, all three codes predict the leeward Cp well but larger deviations are seen on the windward 
side. On the windward side, OVERFLOW underpredicts for a = 12° underpredict for a = 20°. 



Figure 19. Flow visualizations of run 307, a 0°, 0 0°, quad nozzle, Ct 2. 
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OVERFLOW FUN3D DPLR TEST OVERFLOW FUN3D DPLR TEST 




Figure 20. Flow visualizations of run 307, ex. 12°, <fi 0°, quad nozzle, Ct 2. 




Figure 21. Flow visualizations of run 307, ex. 20°, < p 0°, quad nozzle, Ct 2. 
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(a) Model surface Cp 



(b) Model C A , total, C N: and C m 

Figure 22. Comparison of averaged CFD and test surface pressure coefficient as well as code-to-code comparisons of 
C a, total, Cn , and Cm as a function of time. Run 307, quad nozzle, Ct 2, <fi 0°. 
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V.F. Run 311, Quad Nozzle, Cp 2, </> 180' 


The same steadier behavior seen in run 307 at a = 0° was seen in run 311 at a = 0°, as expected. This 
behavior was properly simulated by DPLR and FUN3D. OVERFLOW predicted an unsteady bow shock 
behavior that oscillated between near correct and very large standoff distances. Even with the differences 
in behavior, the Cp on the model face was similar to that seen for run 307 at a = 0°. On the side of 
the model OVERFLOW predicted lower pressure than DPLR and FUN3D, but were mostly still inside the 
tunnel uncertainty. 

For a = 12°, a steadier mode was observed where the windward bow shock was somewhat stable, and the 
leeward bow shock showed an unsteady oscillatory shedding. This behavior was predicted most accurately 
by DPLR. Toward the end of the test Schlieren video, this mode shifted to more like the behavior seen in 
the OVERFLOW results. 

An even steadier mode was observed in the a = 20° than the a = 12° case. FUN3D and OVERFLOW 
capture this mode while the DPLR results were unsteady. All codes predict the surface Cp well except for on 
the windward face between the peripheral nozzle and shoulder, and just aft of the shoulder on the windward 
side. 



Figure 23. Flow visualizations of run 311, ct 0°, <fi 180°, quad nozzle, Ct 2. 
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Figure 24. Flow visualizations of run 311, at. 12°, <fi 180°, quad nozzle, Ct 2. 


Time 1 Time 2 Time 3 Time 4 




Figure 25. Flow visualizations of run 311, ex. 20°, <fi 180°, quad nozzle, Ct 2. 
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(a) Model surface Cp 



(b) Model C a, total, C N , and C m 

Figure 26. Comparison of averaged CFD and test surface pressure coefficient as well as code-to-code comparisons of 
C a, total, C]\r, and Cm as a function of time. Run 311, quad nozzle, Ct 2, </> 180°. 
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VI. Discussion of Results 


One of the largest differences in the code-to-code comparison was the level of unsteadiness predicted by 
each solver. DPLR, FUN3D, and OVERFLOW all used the SST turbulence model, but the implementation 
was slightly different for all. DPLR used the vorticity based production term, while OVERFLOW used the 
strain-based production term with the realizability constraint. OVERFLOW used DES with SST as the 
submodel for the single and triple nozzle cases, and RANS SST for the zero and quad nozzle cases. FUN3D 
used DES with SA as the submodel. Each of these turbulence models generates different levels of eddy 
viscosity, which may add dissipation to the solutions making them more steady. Generally speaking, the 
RANS models were more steady than the DES and the voriticity based SST models were more steady than 
the strain based with realizability constraint, although grid resolution and code implementation may also 
significantly contribute to those trends. 

Grid resolution is also another difference among the flow solvers. Even though refinement for all grids 
is focused on the plume region, differences in spacing and topology do exist (see reference [9]). The level of 
unsteadiness is influenced by the level of dissipation, which can be increased by the coarseness of the grid. 

The ability to properly capture the unsteady effects had varying effect on matching test averaged surface 
Cp. Solving the simulations time-accurately took large computing resources to obtain an average. It was 
shown that in some circumstances a steady state solution compared as well as a time-accurate one (for 
example, run 262 a = 16° or run 311 a = 12°). Other cases, it seemed to matter that unsteady effects were 
captured (for example run 165 a = 0°, a = 12° or run 262 a = 12°). No specific trend was established 
meaning a case by case inspection would be required. Since that is the case, to be able to predict an SRP 
flowfield where test or flight data is not available, a time-accurate approach would be required since it would 
not be possible to know if a steady state solution would be sufficient. 

Obtaining a proper average from such an unsteady flowfield was an issue with the CFD results. A small 
time step was required to resolve the high frequency unsteadiness, and yet a large amount of time is needed 
to define an average. The periodic unsteady behavior seen in run 165 at a = 0° was more conducive to 
calculating a proper average, but for other cases, the unsteadiness was random, and no periodicity was noted. 
It may have been possible for the CFD to resolve a certain section of the unsteady behavior and obtain a 
converged average, but still not cover the entire 2.5 seconds of the test data acquisition window, and not 
obtain the same average as the pressure taps. 

For the wind tunnel, the pressure taps (not including the high frequency gauges) took 25 or 75 readings 
over 2.5 seconds, giving rates of 10 or 30 readings per second. For the CFD results, total run times were 
under 0.05 seconds with thousands of data points for averaging. The wind tunnel sample interval (0.1 or 
0.033 seconds) is on the order or greater than the total CFD run times. This means flowfield variations with 
time scales longer than 0.1 seconds would be averaged out in the wind tunnel data but totally missed in the 
CFD. Also, flowfield variations less than the tunnel sample interval were captured and averaged by the CFD 
but totally missed by the wind tunnel data. 

With the inherent unsteadiness of the flowfield, vehicle stability becomes a concern. It is important to 
realize that the measured unsteadiness in the surface tap readings have a small contribution to the total axial 
force. Figure 27 shows a bar chart depicting the contributions of aerodynamics and thrust to the overall 
axial force coefficient. In most cases, the surface aerodynamic forces contribute less than 5%. This should 
result in the surface pressures creating an even smaller percentage of the total axial force. Also, the required 
flight thrust coefficients will be much higher (~10s) than those achieved in the test. Thrust coefficients this 
large were not obtained in this test due to tunnel interference. A couple runs were made for the triple nozzle 
configuration at a = 0° with Cp values as high as 6. At this higher thrust coefficient, the flow was more 
steady as seen in Figure 28 and the contribution of aerodynamic effects to the total axial force was much 
less than for smaller thrust coefficients. Future plans of the project is to take the same 5 inch diameter 
model into the Ames 9’x7’ Unitary Plan Wind Tunnel. With the larger test section, tunnel wall effects will 
be minimized and larger thrust coefficients will be possible. 
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Figure 27. Contributions of aerodynamic surface forces and thrust to the total axial force coefficient. 



Figure 28. Wind tunnel Schlieren (left) and constructed Schlieren of OVERFLOW simulation (right) of the triple 
nozzle configuration at Mach 3.5 and Ct = 6. This higher thrust coefficient demonstrated steadier behavior. 


VII. Conclusions 

Computational fluid dynamics has taken a strong step towards validation for supersonic retropropulsion. 
Using four solvers, DPLR, FUN3D, and OVERFLOW, comparisons were made to test 1853 from the Lang- 
ley supersonic Unitary Plan Wind Tunnel. Code-to-code and code-to-test comparisons are encouraging and 
possible error sources have been defined. An inherent unsteadiness was noted at the tested thrust coeffi- 
cients, but it was noted that the flowfield is much steadier at higher thrust coefficients for the triple nozzle 
configuration. 
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